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Abstract 

We  introduce  simplex  free  adaptive  tree  numerical  methods  for  solv¬ 
ing  static  and  time  dependent  Hamilton- Jacobi  equations  arising  in  level 
set  problems  in  arbitrary  dimension.  The  data  structure  upon  which  our 
method  is  built  is  a  generalized  n- dimensional  binary  tree,  but  it  does 
not  require  the  complicated  splitting  of  cubes  into  simplices  (aka  gen¬ 
eralized  n-dimensional  triangles  or  hypertetrahedrons)  that  current  tree 
based  methods  require.  It  has  enough  simplicity  that  minor  variants  of 
standard  numerical  Hamiltonians  developed  for  uniform  grids  can  be  ap¬ 
plied,  yielding  consistent,  monotone,  convergent  schemes.  Combined  with 
the  fast  sweeping  strategy,  the  resulting  tree  based  methods  are  highly 
efficient  and  accurate.  Thus,  without  changing  more  than  a  few  lines  of 
code  when  changing  dimension,  we  have  obtained  results  for  calculations 
in  up  to  n  =  7  dimensions. 


1  Introduction 

In  this  paper  we  present  a  simplex  free  adaptive  tree  numerical  method  for 
solving  static  and  time  dependent  Hamilton- Jacobi  partial  differential  equa¬ 
tions  (H-J  PDEs)  arising  in  level  set  problems  in  arbitrary  dimension.  The 
method’s  adaptivity  increases  resolution  near  the  interface  being  studied,  and 
simplifies  previous  successful  tree  based  implementations,  allowing  for  its  exten¬ 
sion  to  arbitrary  dimension  without  an  increase  in  the  complexity  of  function 
reconstruction,  which  is  a  necessary  part  of  finding  spatial  derivatives  needed 
in  solving  the  PDEs. 
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Applications  in  higher  dimensions  requiring  adaptive  meshes  to  resolve  fine 
details  arise  in  numerous  fields.  In  [23], [15], [9], [29]  multi-valued  solutions  to  H- 
J  equations  were  found  by  replacing  a  single  valued  solution  with  the  level  set 
(or  intersection  of  level  sets)  of  a  higher  dimensional  function.  This  idea  was 
also  used  in  [2], [6], [8]  to  study  interfaces  with  codimension  >1.  In  [32], [33]  a 
level  set  formulation  was  used  to  solve  problems  arising  in  mathematical  finance, 
where  high  dimensional  issues  are  routinely  encountered.  Even  for  codimension- 
1  problems  in  3d  there  is  still  a  desire  to  find  implementable  adaptive  methods 
to  resolve  fine  details,  such  as  in  the  segmentation  of  the  human  brain,  or  other 
applications  involving  highly  curved  surfaces  such  as  Wulff  crystals.  See  [24], [31] 
for  a  wide  range  of  physical  problems  to  which  level  set  methods  are  applied. 

Since  the  introduction  of  level  set  methods  for  interface  tracking  [22]  there 
has  been  work  done  in  an  attempt  to  reduce  the  component  of  the  computa¬ 
tional  portion  of  the  method  subject  to  the  most  criticism:  the  necessity  of 
extra  dimensions.  Within  a  few  years  following  [22]  narrow  band  methods  were 
proposed  that  reduced  the  computational  complexity  by  resolving  the  level  set 
function  only  near  the  interface  being  tracked  [1] ,  [38] ,  [26] .  These  methods  were 
able  to  use  the  well  established,  convergent,  finite  difference  schemes  available 
to  uniform  grids. 

However,  these  narrow  band  methods  did  not  reduce  the  storage  require¬ 
ments,  limiting  them  to  the  same  resolutions  which  uniform  grids  were  re¬ 
stricted.  Following  this,  tree  based  methods  were  introduced,  allowing  for  adap¬ 
tivity  of  the  mesh  near  the  interface,  while  not  sacrificing  too  much  complexity 
[34], [21], [11], [19].  The  tree  data  structure  used  in  these  methods  was  well  under¬ 
stood  by  the  computer  science  community,  and  thus  data  storage  and  retrieval 
were  able  to  be  carried  out  in  an  efficient  manner.  However,  the  nonuniformity 
of  the  mesh  required  new  schemes  to  be  developed  for  the  PDEs  to  be  solved. 
In  some  cases  semi-Lagrangian  CIR  [10]  schemes  were  used  for  time  dependent 
level  set  equations.  These  schemes  have  some  drawbacks,  though.  Firstly,  they 
are  only  provably  convergent  for  hyperbolic  problems,  and  many  level  set  PDEs 
involve  mean  curvature  or  are  otherwise  parabolic  in  nature.  Secondly,  they  re¬ 
quire  a  backtracking  along  characteristics  and  an  interpolation  at  an  arbitrary 
point  within  the  domain.  This  interpolation  is  a  delicate  process  that  requires 
the  division  of  the  domain  into  simplices,  which  can  become  complicated  in 
higher  dimensions  [21],  In  [19]  CIR  was  used  for  advection  of  values  stored 
at  cell  corners,  and  cell  centered  data  was  stored  for  the  pressure  equation  in 
Navier-Stokes,  where  a  one  point  (constant  within  each  cell)  interpolation  tech¬ 
nique  was  used  to  avoid  apparent  complexities,  and  to  preserve  the  symmetry 
of  the  discretization.  In  addition,  for  the  eikonal  equation  used  to  maintain 
the  signed  distance  property  of  the  level  set  function,  [19]  used  some  special 
treatment  at  T-junctions  in  the  context  of  the  fast  marching  spirit  [37] . 

There  have  been  other  local  level  set  methods  [20], [35], [5], [14], [4]  proposed 
which  range  from  variants  of  AMR  to  using  tubes  of  uniformly  spaced  grid 
points  near  the  interface.  Some  of  the  methods  approach  the  complexity  of  [26] , 
eliminating  the  need  to  store  the  unused  grid  points  away  from  the  interface  of 
interest.  They  also  allow  for  the  standard  finite  difference  schemes  to  be  used 
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as  the  grid  is  uniform  near  the  interface.  However,  with  these  gains  comes  addi¬ 
tional  complexity  in  implementation,  and  it  should  be  noted  that  the  successive 
improvements  and  acceptance  of  the  tree  based  methods  in  various  applications 
are  testaments  to  their  facility  and  usefulness. 

In  this  paper  we  introduce  a  tree  based  method  that  retains  the  advantages 
of  the  previous  tree  based  algorithms,  such  as  having  a  well  studied  and  un¬ 
derstood  data  structure,  while  avoiding  the  drawbacks  of  having  inconsistent 
schemes  requiring  n-dimensional  simplices  and  interpolation.  Thus  we  are  able 
to  use  the  standard  numerical  Hamiltonians  derived  for  uniform  grids  (modified 
slightly)  which  result  in  consistent,  monotone,  convergent  numerical  methods. 
Combined  with  the  fast  sweeping  strategy  [40], [36], [16], [28],  the  resulting  tree 
based  methods  are  highly  efficient  and  accurate. 

The  paper  consists  of  a  brief  overview  of  the  tree  data  structure,  followed  by 
a  discussion  of  the  numerical  schemes  for  static  H-J  equations,  and  then  time 
dependent  H-J  equations.  Finally,  numerical  results  are  given  for  codimension- 1, 
codinrension-2,  and  codimension-n  problems. 


2  Tree  Data  Structure 

In  this  section  we  describe  the  tree  data  structure  used.  We  use  a  generalized 
binary  tree  (e.g.  quadtree  in  2d,  octree  in  3d,  etc.)  data  structure,  details 
of  which  can  be  found  in  numerous  computer  science  texts  [18], [30], [13].  We 
describe  the  portions  of  the  implementation  that  are  specific  to  our  problem  of 
solving  a  PDE  in  a  bounded  spatial  domain. 

We  assume  a  computational  domain,  f l  =  [0,  l]n.  At  the  kth  level  of  the  tree, 
each  node,  c,  represents  a  hypercube  cell  with  sides  of  length  dxc  =  l/2fc,  and 
center  xc.  We  assume  that  the  level  set  function  value,  <f>,  is  stored  at  the  centers 
of  mass  of  the  nodes  at  the  finest  level  of  the  tree,  also  known  as  the  leaves. 
When  refinement  of  a  cell  is  done,  the  cell  is  split  into  2"  subcells  with  side 
lengths  l/2fc+1.  We  do  not  allow  any  cells  with  side  length  ratio  >  2  or  <  0.5  to 
be  neighbors.  This  final  restriction  can  be  obtained  by  following  the  criterion  of 
refining  any  cell  whose  distance  to  the  interface,  T,  is  less  than  a  constant  times 
its  edge  length  [34].  In  practice,  if  we  are  using  a  single  level  set  function  <j)  e.g. 
for  codimension- 1  problems,  if  <j>  is  a  signed  distance  function  then  we  can  set 
p  >  (l  +  y/n/2)  and  refine  if  |^(a:c)|  <  pdxc.  For  problems  where  the  intersection 
of  multiple  level  set  functions,  {<j>j},  represents  T,  where  the  level  sets  of  the  4>j 
are  mutually  orthogonal  and  each  <pj  is  a  distance  function  measured  along  the 
level  sets  of  the  other  then  we  can  compare  \\(f>\\l2  to  pdxc. 


3  Static  H-J  Equations 

Here  we  introduce  a  fast  sweeping  implementation  for  solving  certain  static 
Hamilton- Jacobi  equations  such  as  the  eikonal  equation  |V^|  =  /  or  in  general 
H(V<j>)  =  f  [36], [39], [40], [17].  These  type  of  equations  are  commonly  found 
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when  reinitializing  the  level  set  function  to  be  a  signed  distance  function  during 
a  dynamic  evolution,  or  for  other  weighted  distance  calculations  that  arise  in 
numerous  physical  problems. 

In  order  to  avoid  n-triangulations  that  can  lead  to  very  complicated  local 
Hamiltonian  solvers  [28] ,  we  follow  the  same  ideology  that  many  other  adaptive 
and  tree  based  methods  follow:  study  a  small  number  of  local  node  config¬ 
urations  of  the  grid,  and  then  appropriately  scale  them  so  that  the  various 
operations  of  interpolation,  refinement,  etc.  can  be  applied  in  the  same  way 
anywhere  in  the  domain. 


3.1  Local  Numerical  Hamiltonian  Solver 

The  first  part  of  the  fast  sweeping  method  is  the  local  numerical  Hamiltonian 
solver.  We  will  design  monotone,  consistent  solvers  that  do  not  require  n- 
triangulations. 


3.1.1  Upwind  Hamiltonians 

The  upwind  Hamiltonians  approximating  the  eikonal  equation  |Vit|  =  /  with 
boundary  data  given  on  T  are  based  on  using  the  Godunov  Hamiltonian  (GH): 


H(D*}u(y),  Dx+'u(y), . . . ,  Dx_"u(y),  Dx+"u{y))  = 


n 


ma x{(U?«(y))+,(U^W(y))-}2 


(1) 


where  y  is  a  grid  point  where  we  wish  to  update  the  numerical  solution  u1  or 
the  Osher-Sethian  Hamiltonian  (OSH): 


H(DxJu(y),Dx+'u(y), . . . ,  Dx_"u(y),  Dx+”u(y))  = 


,|E[(^-*«(»))+]2  +  [W«(l/))-]2-  (2) 

\  i=i 


Godunov  numerical  Hamiltonians  were  evaluated  in  [3], [25].  For  nonlinear  HJ 
equations  whose  Hamiltonians  differ  significantly  from  that  of  the  eikonal  equa¬ 
tion  the  resulting  expressions  become  quite  complicated,  involving  many  “if” 
statements. 

When  the  grid  is  uniform  along  the  ith  axis,  the  standard  2  point  finite 
difference  can  be  used  for  e.g.  DxJu{y ),  by  taking  u(y)-u(y~Sei) ;  where  8  is  the 
local  spacing  between  nodes  in  the  ith  direction. 

Admittedly,  the  grid  is  not  uniform  everywhere,  so  when  we  try  to  compute 
quantities  such  as  D^ufy)  we  will  not  have  a  uniform  definition  throughout  all 
of  the  domain.  However,  because  the  grid  refinement  is  done  predictably  (by 
this  we  mean  that  there  are  only  a  small  number  of  local  configurations  of  the 
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grid,  up  to  scaling),  we  can  quickly  find  the  value  at  an  offset  point  in  the  —  ej 
direction,  u(y  —  Se/),  needed  in  D°^u(y). 

We  assume  that  the  tree  is  constructed  so  that  in  each  cell  S,  us  is  defined 
at  the  center,  ys,  of  S,  and  each  cell  is  an  ?i-cube  with  equal  side  lengths  given 
by  5s  ■  This  uniform  restriction  can  be  relaxed,  but  for  expositional  purposes  it 
will  not  be.  Also,  we  note  the  restriction  that  the  ratio  of  the  5  of  neighboring 
cells  is  either  2,  1  or  1/2. 

In  any  dimension,  n,  there  are  only  3  possible  local  configurations  that  need 
examining,  when  attempting  to  find  D^u(y): 

1.  The  cell  B  that  is  adjacent  to  cell  A  9  y  in  the  —  e*  direction  is  exactly 
the  same  size  as  A,  thus  u(y  —  Sei)  =  Ub,  5  =  5a  =  5b-  This  is  standard 
2  point  finite  differencing. 

2.  The  cell  B  that  is  adjacent  to  cell  A  9  y  in  the  —ei  direction  is  smaller 

than  A ,  i.e.  5b  =  5a/ 2.  In  this  case  we  have  a  situation  illustrated  in 
Figure  1  in  2d.  In  n  dimensions  we  take  u(y  —  <5ej)  =  {  the  average  of 
the  2n  adjacent  neighboring  cells  whose  faces  with  outward  normal  e*  are 
touching  the  face  of  A  that  has  outward  normal  —  }.  As  the  centers  of  all 

these  points  are  coplanar,  this  is  just  the  linearly  interpolated  value  at  the 
point  P  that  is  the  intersection  of  the  plane  containing  these  neighboring 
cell  centers  and  the  ray  given  in  parametric  form  by  ?’(r)  =  y  —  rei ,  r  >  0. 

3.  The  cell  B  that  is  adjacent  to  cell  A  9  y  in  the  — e*  direction  is  larger  than 
A ,  i.e.  5b  =  “25a-  In  this  case  ys  does  not  lie  along  r(r).  However,  because 
of  the  structure  of  the  grid  points  in  the  tree,  there  is  a  neighboring  cell,  C, 
of  A  such  that  yWyc  intersects  r(r).  This  cell  C  is  the  diagonal  neighbor 
of  A  in  the  direction 


(sgn(yA, i  -  Vb, i),  •  •  • ,  sgn(yA,i- i  -  2/s.i-i), 

-  sgn(yA,i  - 

sgn(yA,i+ 1  -  VB,i+ 1),  •  •  • ,  sgn(yA,n  ~  VB,n )), 


where  sgn  is  the  signum  function. 

In  this  case  there  are  2  possibilities  for  C:  case  1.  C  is  either  the  same 
size  as  A;  case  2.  C  is  the  same  size  as  B.  See  Figures  2,  3  for  diagrams 
of  these  cases  in  2d  and  3d. 

The  interpolated  value  at  the  intersection  point  P  is 


u(P) 


1 

\yWyc\ 


(u(vb)  \vbP\  +u{yc)  \yc-P\ ) 


u{yB)wB  +  u{yc)wC- 

(3) 


Also,  we  have  5  =  wB  \yB,i  -  yA,i\  +  wc  \yc,i  ~  VA,i\-  In  case  1  we  find 
Wb  =  2/3,  wc  =  1/3,  and  in  case  2  we  find  Wb  =  3/4,  wc  =  1/4.  These 
are  the  weights  for  any  dimension  n,  which  is  very  appealing  in  that  we 
do  not  have  to  resort  to  complicated  n-triangulations. 
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P 

»  -  -  _  - 

/ 

B 

C 

Figure  1:  Case  where  neighboring  cells  are  smaller  than  A. 
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C 

°  r 

B 

\ 
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Figure  2:  Cases  where  neighboring  cell  B  is  larger  than  A  in  2d.  Left:  case  1, 
right:  case  2. 


We  define  <5j;_  =  S  in  this  particular  case  because  ys,i  <  UA,i-  In  the  case 
where  ys,i  >  VA,i ,  5»,+  is  defined  as  the  distance  between  P  and  A. 

Once  one  has  found  all  the  D±u(y ),  \/i  one  can  solve  the  quadratic  equation 
arising  from  the  local  numerical  Hamiltonian  for  u(y),  and  update  the  solution 
with  this  found  value.  Because  the  <5;,±  in  each  particular  D±u(y)  could  be 
different,  the  Godunov  solver  introduced  in  [40]  with  its  simple  min  and  if 
statements  is  not  applicable.  However,  the  procedure  for  finding  the  correct 
solution  of 

'(x  -  ai)+ 
hi 

can  be  used.  We  present  the  case  for  OSH,  as  GH  is  more  complicated  (but 
feasible) . 

1.  Let  the  a,j ,  j  =  1  :  2 n  be  the  points  from  the  finite  differences 
{dj}  =  {u(yA  ±  eAi±)},  i  =  l:n, 

ordered  from  least  to  greatest,  and  the  hj  be  the  corresponding  offset 
distances  from  yA  (the  hj  will  not  be  in  any  particular  order).  We  set 

«2n+l  =  oo. 


- h 


(x  -  am)+ 


=  f 


(4) 
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Figure  3:  Cases  where  neighboring  cell  B  is  larger  than  A  in  3d.  Left:  case  1, 
right:  case  2. 


2.  Set  m  =  1; 

3.  Solve  J^jLi [^"^  ]2  —  /  t°  get  a  solution  x. 

4.  Check  to  see  if  x  <  am+ 1.  If  so,  then  we  are  done  and  we  set  u(y)  =  x.  If 
not,  then  we  set  m  — >  m  +  1,  and  go  to  step  3,  unless  rn  =  2 n,  then  we 
are  done. 

Note:  the  implementation  of  this  solution  algorithm  is  independent  of  n, 
except  for  the  number  of  terms  in  the  sum. 

The  algorithm  for  GH  is  more  complicated  as  GH  includes  max  functions 
that  OSH  does  not.  The  only  drawback  of  OSH  is  its  slightly  larger  error  near 
sonic  shocks,  but  if  we  refine  the  grid  such  that  it  is  locally  uniform  at  sonic 
shocks,  then  GH  could  be  used  there  (and  anywhere  else  on  the  grid  that  is 
locally  uniform),  as  it  can  be  solved  by  the  method  presented  in  [40]. 

Note  that  monotonicity  is  satisfied  because  of  the  positive  weights,  w,  mul¬ 
tiplying  the  u(z),z  ^  y  in  each  finite  difference.  Also,  because  of  the  linear 
interpolations  used,  the  scheme  is  consistent.  Thus,  the  scheme  is  convergent. 

3.1.2  Lax-Friedrichs  Hamiltonian 

Here  we  present  a  Lax-Friedrichs  Hamiltonian  (LFH)  which  does  not  require 
nonlinear  inversions  when  it  is  being  solved.  This  extends  its  applicability  to  a 
wide  range  problems  including  those  with  nonconvex  Hamiltonians.  This  is  a 
generalization  of  the  Hamiltonian  presented  in  [16]. 

The  numerical  Hamiltonian  for  H(Vu)  =  f  with  boundary  data  given  on  F 
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is  as  follows: 


H(DxJu(y),D^u(y), D*_"  u(y) ,  D?  u(y))  = 

H{wf Dx_fu(y)  +  wf  Dx^u{y), . . . ,  w~Dx_nu(y)  +  w+Dx+nu(y)) 

n 

-  ai(D+u(y)  ~  D~u(y)),  (5) 


where 


wf  = 


Ji,+ 


w,  = 


K- 


<5?:,+  +  &i,—  1  <5?:,+  +  di,- 


(6) 


Note  that  +  wi  =1  so  the  scheme  is  consistent. 

To  determine  the  size  of  cr,;  we  note  that  monotonicity  requires  that 

dH/dp+  <  0,  d H/dp~  >  0. 

This  leads  to  the  requirement 

o’,:  >  \HPi\max(wf,w^). 

Within  the  fast  sweeping  framework,  in  order  to  advance  the  solution  a  single 
iteration,  one  writes  (5)  in  the  form 

H  =  L(u(Vl  \  y))  +  cu(y), 

and  then  the  advancement  can  be  written  as 


un+\y)  = 


f(y)-L(un(n\y)) 


The  weights  w  are  composed  precisely  so  that  the  coefficient  c  can  be  calculated 
in  a  linear  way  purely  from  the  artificial  diffusion  term. 

Implementation  of  this  LFH  is  simpler  than  GH  or  OSH  because  it  does  not 
require  the  inversion  of  H .  Thus  is  does  not  require  any  if  statements  and  is 
thus  faster  per  iteration.  However,  it  does  require  more  iterations  to  converge. 


3.2  Sweeping  Directions 

The  second  part  of  the  fast  sweeping  method  is  to  determine  the  directions 
of  the  sweep.  In  [28]  a  strategy  using  reference  points  was  introduced,  with 
an  initial  sequential  ordering  of  the  nodes  with  respect  to  their  distances  from 
these  points.  This  could  work  for  our  method,  but  the  built  in  structure  of  the 
tree  allows  for  another  method  of  sweeping. 

The  tree  specific  sweeping  method  is  as  follows: 

1.  Each  ordering  of  sweeping  is  defined  by  an  ordering  of  the  vertices  of  the 
n-cube.  In  n  dimensions  there  are  2n  possible  sweeps  that  are  found  by 


A 

B 

C 

D 

Figure  4:  Sample  child  ordering  in  2d.  The  outer  perimeter  is  the  boundary  of 
the  parent  cell.  Each  lettered  interior  square  corresponds  to  a  child  pointer  to 
one  of  the  4  smaller  squares. 


taking  all  combinations  of  sweeping  from  low  to  high,  or  high  to  low  in 
each  dimension.  So  in  2d  if  the  vertices  of  a  square  are  as  in  Figure  4  the 
possible  orderings  are 


{A,B,C,D}, 

{D,C,H,H}, 

{C,D,A,  B}, 

{B,A,D,C}. 

2.  For  each  of  the  orderings  in  step  1,  call  a  preorder  traversal  [18]  of  the  grid 
starting  at  the  root  node,  based  on  a  particular  ordering  of  the  children 
given  from  the  previous  step.  When  a  leaf  is  reached,  update  the  solution 
using  the  local  Hamiltonian  solver. 

This  means  that  if  we  choose  the  child  ordering  {A,  B,C,  D}  then  we  call 
a  preorder  traversal  with  node  A  as  the  starting  node,  followed  by  a  preorder 
traversal  with  node  B  as  the  starting  node,  etc.  Once  the  traversal  has  gotten 
to  the  leaf  of  the  tree  (i.e.  a  node  with  no  children)  we  update  u.  This  recursive 
type  of  tree  visitation  is  standard  and  can  be  found  in  any  thorough  book  on 
computer  algorithms  and  binary  trees  [18], [30]. 

Figure  5  shows  a  sample  ordering  of  the  nodes  when  all  nodes  are  at  a 
uniform  depth  in  the  tree.  In  this  particular  case  the  sweeping  algorithm  will 
give  exactly  the  same  result  as  the  standard  sweep  ordering  [40]  for  a  uniform 
grid  of  this  size  with  a  fixed  node  at  the  origin  with  u( 0,  0)  =  0.  This  is  because 
for  each  node  visited  in  this  sweep,  all  nodes  to  the  south  and  west  of  it  have 
already  been  updated  in  the  sweep. 

However,  when  the  grid  is  not  uniform  we  have  found  that  extra  sweeps  are 
needed  as  D^u^)  may  depend  on  nodes  that  are  farther  north  or  east  than 
y ,  as  is  the  case  when  the  stencil  choices  in  Figure  2  would  be  used.  We  make 
some  comments  on  this.  Firstly,  since  the  sweeping  strategy  accesses  all  nodes 
systematically,  the  tree  based  methods  will  converge  eventually,  no  matter  how 
many  sweeps  it  needs.  Secondly,  since  there  are  many  more  possible  information 
flowing  directions  on  a  non-uniform  grid  as  demonstrated  in  [28],  and  the  tree 
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Figure  5:  Nodes  that  are  leaves  of  a  tree  on  a  uniform  16  x  16  grid,  and  the 
path  connecting  sequential  nodes  in  the  sweep  from  bottom  left  to  top  right. 


specific  sweeping  method  designed  here  will  only  allow  a  finite  number  of  such 
information  flowing  directions  to  be  treated  simultaneously,  it  is  reasonable  for 
the  tree  specific  methods  to  use  extra  sweeps  to  converge.  Thirdly,  it  is  possible 
to  design  optimal  tree  specific  sweeping  methods  by  reordering  the  nodes  so 
that  all  the  possible  information  flowing  directions  can  be  covered  with  a  finite 
number  of  tree-based  orderings. 

3.3  Solution  Procedure 

The  complete  solution  procedure  is  as  follows: 

1.  Assume  we  are  given  a  grid  upon  which  the  solution  <j>  will  be  solved.  Find 
all  points  within  a  small  0(dx )  distance  of  F  and  find  an  approximate 
solution  (f>n  at  these  points.  This  can  be  done  by  interpolating  the  known 
boundary  F.  Set  all  other  points  to  a  large  value,  e.g.  <pn  =  1010,  which 
is  larger  than  the  exact  solution  on  the  grid. 

2.  Sweep  through  the  domain  {xj}  using  a  preorder  traversal  specified  by 
one  of  the  2n  orderings  of  child  pointers,  solving  for  (j>*{xj)  at  each  grid 
point  using  Gauss-Siedel  iteration  by  inverting  GH,  OSH  or  LFH.  Take 

(, bn+ 1  =  min{4>*{xj),4>n{xj)). 

3.  Check  if  \\(/>n  —  </>"+1||oo  <  G  where  e  is  a  fixed  tolerance  to  indicate  con¬ 
vergence  has  been  reached.  If  convergence  has  not  been  reached  go  to  step 

2. 
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4  Time  Dependent  Level  Set  Equations 

This  section  concerns  solving  time  dependent  H-J  equations  as  well  as  higher 
order  parabolic  equations  such  as  mean  curvature  motion  that  arise  frequently 
in  level  set  problems.  We  will  introduce  the  way  in  which:  1.  the  numerical 
Hamiltonian  is  constructed;  2.  new  values  are  chosen  on  the  grid  in  a  monotone 
way  after  a  refinement/coarsening  has  taken  place. 


4.1  Numerical  Hamiltonian 


Examining  the  constructions  for  the  functions  at  the  points  y  ±  <5e,;  from  the 
section  on  static  H-J  equations  we  can  see  that  2  point  upwind  differences  can 
be  calculated  easily  for  any  fine  grid  point  y.  For  example 


DxJu(y) 


u(y)  -  ( u(yB)wB  +  u(yc)wc) 

6 


(7) 


using  the  notation  from  (3).  Then  standard  monotone  numerical  Hamiltonians 
such  as  GH,  OSH  and  LFH  can  be  used  with  Runge-Kutta  (R-K)  timestepping 
with  the  CFL  condition  determined  by  the  finest  cell  size  being  used.  For 
higher  order  WENO  type  reconstructions  a  bit  more  work  would  be  involved, 
but  they  are  certainly  possible  to  construct  and  we  would  still  avoid  the  need 
to  use  high  dimensional  triangulations.  Central  differences  and  higher  order 
derivatives  can  be  calculated  by  using  the  weightings,  w,  that  were  introduced 
in  section  3.1.2  for  Lax-Friedrichs  Hamiltonians  after  calculating  the  first  order 
upwind  differences  at  the  grid  points. 


4.2  Interpolation  After  Refinement/Coarsening 

We  would  like  the  interpolation  processes  to  be  consistent  and  monotone  as  well 
so  that  our  entire  scheme  including  adaptation  is  stable.  Also,  we  would  like  to 
avoid  having  to  work  with  any  high  dimensional  simplices.  It  turns  out  that  as 
was  the  case  with  the  upwind  differencing,  we  have  only  a  few  different  cases 
that  can  be  applied  to  all  dimensions  in  the  same  way. 


4.2.1  Coarsening 

After  a  coarsening  the  value  at  the  center  of  the  node  of  size  l/2k  that  just 
became  a  leaf  is  set  to  be  the  average  of  the  2"  function  values  at  the  nodes 
with  side  lengths  l/2fe+1  that  were  its  children.  Coarsening  is  done  if  all  the 
siblings  of  a  cell  c  are  at  the  same  level  as  c  and  have  not  been  marked  for 
refinement,  and  if  the  resulting  averaged  coarsened  value,  0c°arsened;  does  not 
violate  the  refinement  condition  (i.e.  we  do  not  have  |  <^coarsened  |  <  2 pdxc). 

4.2.2  Refinement 

Examining  Figure  6  we  can  see  that  after  refinement  we  can  use  2  point  linear 
interpolation  to  find  the  new  value  at  point  A  with  weights:  in  case  1,  wb  = 
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B 

\  A 

B 

\  A 

Figure  6:  Possible  refinement  cases  in  2d.  A:  Node  where  interpolation  is  eval¬ 
uated.  B:  Node  at  center  of  coarse  cell  that  was  divided.  C:  Diagonal  neighbor 
used  in  interpolation.  Left:  Case  1.  Right:  Case  2. 


2/3,  wc  =  1/3;  and  in  case  2,  wb  =  3/4,  wc  =  1/4.  These  are  the  weights 
for  any  dimension  where  the  diagonal  neighbor,  C,  is  the  adjacent  cell  in  the 
direction 


(sgn(yA, l  -  Vb, i),  •  •  • ,  sgn(yA,n  ~  VB,n))- 

Both  the  coarsening  and  refinement  interpolations  are  monotone  and  con¬ 
sistent  for  linear  functions  </>. 

4.3  Solution  Procedure 

The  complete  solution  procedure  for  the  time  dependent  problem  is  as  follows: 

1.  Assume  we  are  given  a  grid  upon  which  the  solution  (/>  will  be  solved,  and 
we  are  given  initial  conditions  cj)n .  If  />"  is  not  a  signed  distance  function 
then  use  the  fast  sweeping  method  with  OSH  or  GH  to  reinitialize  the 
solution. 

2.  For  each  time  step,  for  every  grid  point,  advance  the  time  dependent  H-J 
equation  forward  one  time  step. 

3.  Reinitialize  the  solution  using  fast  sweeping. 

4.  Refine  the  solution  where  necessary. 
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5.  Coarsen  the  solution  where  necessary. 

6.  Go  to  step  2. 

It  should  be  noted  that  the  reinitialization/refinement/coarsening  proce¬ 
dures  do  not  need  to  be  carried  out  every  time  step,  but  should  be  done  at 
least  once  every  few  time  steps. 

5  Numerical  Results 

In  this  section  numerical  results  are  presented  for  both  static  and  time  dependent 
problems  for  codimension- 1,  codinrension-2,  and  codimension-n  problems.  We 
refrain  from  studying  the  reduced  memory  requirements  that  the  tree  method 
gains  over  uniform  discretizations  as  these  are  well  presented  in  [21].  We  focus 
on  presenting  convergence  rate  estimates  and  show  error  results  in  dimensions 
up  to  n  =  7.  For  the  time  dependent  problems  we  use  forward  Euler  time 
advancement. 

In  Table  1  we  show  errors  and  node  counts  for  a  codimension-n  problem  of 
solving  the  eikonal  equation 


|V0|  =  1,  (8) 

with  a  boundary  point  at  Xb  =  (0.5, 0.5, . . .  ,0.5)  with  value  u{xb)  =  0.  In  this 
table  the  number  of  leaves  represents  the  number  of  points  where  the  function 
value  is  stored,  while  the  number  of  uniform  points  represents  what  this  number 
would  be  if  a  uniform  grid  with  the  finest  dx  listed  was  used.  We  note  that  the 
adaptivity  for  this  problem  is  based  on  knowledge  of  the  distance  to  the  fixed 
node  at  the  center  of  the  domain,  which  is  essentially  knowledge  of  the  solution 
to  (8).  Thus  a  better  adaptivity  routine  needs  to  be  formulated,  which  will  be 
the  subject  of  future  research. 

We  use  the  preorder  traversal  fast  sweeping  method  with  the  OSH  numerical 
Hamiltonian.  The  error  is  measured  on  the  finest  3  levels  of  leaves,  except  in 
7d  where  it  is  measured  on  the  finest  2  levels  of  leaves. 

An  interesting  point  to  note  is  that  when  a  stopping  criterion  is  used  such 
as  stopping  when  the  change  in  the  error  is  <  10~12  we  find  that  the  number 
of  sweeps  needed  to  achieve  this  criterion  does  not  increase  linearly  in  n,  but 
rather  stagnates.  For  example  in  5d  we  need  only  30  sweeps,  in  6d  we  need  35, 
and  in  7d  only  32.  Perhaps  this  is  some  effect  particular  to  our  specific  example, 
but  perhaps  not. 

In  Table  2  we  show  errors  and  convergence  rate  estimates  for  the  time  de¬ 
pendent  H-J  equation 


<k  +  |V0|  =  0,  (9) 

with  initial  condition  given  by  a  hypersphere  of  radius  0.2  centered  at  the  point 
Xb  =  (0.5,  0.5, . . . ,  0.5).  The  normal  motion  moves  the  hypersphere  inwards 
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finest  dx 

L\  error 

Too  error 

N  leaves  of  tree 

N  uniform  points 

o 

1/256 

6.046  E-05 

7.193  E-03 

304 

65,536 

1/256 

1.293  E-05 

1.126  E-02 

2,472 

16,777,216 

□ 

1/256 

3.056  E-06 

1.615  E-02 

26,896 

4.295  E+09 

D 

1/256 

7.209  E-07 

2.027  E-02 

309,536 

1.100  E+12 

H 

1/128 

2.018  E-05 

4.802  E-02 

2,822,464 

4.398  E+12 

□ 

1/16 

1.607  E-01 

2.307  E-01 

9,379,840 

2.684  E+08 

H 

Table  1:  Errors  for  eikonal  equation. 


towards  Xb  with  constant  normal  velocity  =  1.  The  error  is  measured  by  de¬ 
termining  the  volume  of  the  hypersphere  at  the  final  time,  T,  (which  is  taken 
to  be  near  0.1)  and  then  measuring  the  error  in  the  implied  radius  versus  the 
exact  radius  =  0.2  —  T.  This  is  basically  an  L\  error  estimate. 


finest  dx 

error 

rate 

dimension 

1/32 

9.605  E-03 

- 

3 

1/64 

4.376  E-03 

1.134 

3 

1/128 

2.110  E-03 

1.052 

3 

1/256 

1.263  E-03 

0.741 

3 

1/32 

8.201  E-03 

- 

4 

1/64 

4.147  E-03 

0.984 

4 

1/128 

2.138  E-03 

0.956 

4 

1/16 

3.220  E-02 

- 

5 

1/32 

1.112  E-02 

1.535 

5 

1/16 

1.743  E-02 

- 

6 

Table  2:  Convergence  rate  estimate  for  constant  motion  in  normal  direction. 

In  Table  3  we  show  errors  and  convergence  rate  estimates  for  time  dependent 
motion  by  mean  curvature 


4>t  +  k  |V0|  =  0, 


(10) 


where 


k  =  — 


L  i=  1 


E< 

0= 1 


EE< 

*= 1  0= 1 


/|V</>|E 


is  the  scaled  mean  curvature  of  the  set  T  =  {x\(j)(x)}  =  0.  We  initialize  T  as 
a  hypersphere  of  radius  0.2  centered  at  the  point  Xb  =  (0.5, 0.5, . . .  ,0.5).  The 
first  partial  derivatives  of  n  are  found  using  centered  differences  as  suggested 
above  with  the  weights  w  indicated  in  the  section  on  the  LFH,  and  the  second 
derivatives  are  found  using  3  point  centered  stencils. 
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The  mean  curvature  motion  moves  the  hypersphere  inwards  towards  Xb  with 
normal  velocity  =  (n  —  1  )/r(t),  where  r(t)  is  the  radius  of  the  hypersphere  at 
time  t.  The  exact  solution  is  a  hypersphere  with  r(t)  =  \/0.22  —  2 (n  —  l)t, 
centered  at  Xb-  The  final  time  varies  from  dimension  to  dimension  but  we  take 
0(1  / dx giobalmin)  number  of  time  steps  for  the  coarsest  grid  in  each  dimension’s 
convergence  study.  The  error  is  measured  in  the  same  way  as  it  was  for  the 
normal  motion  case. 


finets  dx 

error 

rate 

dimension 

1/32 

1.821  E-03 

- 

2 

1/64 

6.011  E-04 

1.599 

2 

1/128 

1.581  E-04 

1.927 

2 

1/256 

6.457  E-06 

4.614 

2 

1/32 

1.786  E-03 

- 

3 

1/64 

3.217  E-04 

2.473 

3 

1/128 

1.951  E-05 

4.044 

3 

1/16 

2.114  E-02 

- 

4 

1/32 

1.482  E-02 

0.513 

4 

1/64 

7.155  E-03 

1.051 

4 

Table  3:  Convergence  rate  estimate  for  mean  curvature  motion. 


In  Table  4  we  show  error  convergence  rates  for  a  codimension-2  problem  of 
advection  of  a  closed  curve  in  3d  [6] .  We  advance 

(0i)t  +  v  •  V0J- =  0,  V  =  (1,1,1),  (11) 

for  j  =  1,  2,  where  the  initial  condition  is  given  by  a  circle  of  radius  0.2  centered 
at  Xb  =  (0.25,0.25,0.25).  This  circle  is  defined  by  the  intersection  of  the  0  level 
sets  of  2  level  set  functions,  4>i,4>2-  The  final  time  is  T  =  0.3125.  The  exact 
solution  is  a  circle  centered  at  (0.5625,0.5625,0.5625)  with  radius  =  0.2.  For 
this  problem  it  is  necessary  to  set  the  <j>j  to  be  orthogonal  to  each  other  every 
few  time  steps.  This  is  done  in  a  fast  sweeping  way  as  presented  in  [7].  The 
error  is  measured  at  t  =  T  by  sampling  the  circle  that  is  the  exact  solution  with 
500  points,  {y},  and  estimating  the  L\  error  of  \/0i(y)2  +  </>2 (y)2  on  the  circle. 


finest  dx 

error 

rate 

1/16 

4.520  E-02 

- 

1/32 

2.544  E-02 

0.829 

1/64 

1.372  E-02 

0.891 

1/128 

7.261  E-03 

0.918 

Table  4:  Convergence  rate  estimate  for  motion  of  a  closed  curve  in  3d. 
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In  Figure  7  we  show  contour  plots  of  Wulff  shapes  [27]  arising  from  solving 


p2  +  q2  +  r2  +  2|r 


|)  ^1  +  3 ^ 


( p 2  +  q2)3/2  —  (3 p2q  —  q3) 
2  {p2  +  q2)3/2 


=  1. 


(12) 


where  (jp,  q,  r)  =  (<j>x,  <t>y,<t>z).  This  equation  is  equivalent  to  solving 


|V0!  =  (1  +  2|  sin(91|)(l  +  I  sin(1.5(02  +  O.5tt))|)|V0|  =  1, 


where  7  is  known  as  the  surface  tension  in  materials  science,  and  is  a  function 
of  the  spherical  coordinates  (0i,$2)-  For  this  and  the  next  example  we  use 
fast  sweeping  with  LFH.  The  diffusion  terms  satisfy  a  =  (9/4,9/4,27/8).  We 
fix  a  point  in  the  center  of  the  domain  with  the  value  0.  The  grid  is  adapted 
similarly  to  how  it  would  be  for  a  usual  level  set  evolution,  but  with  the  coarsest 
resolution  being  dx  =  1/64,  and  the  finest  resolution,  dx  =  1/1024.  Neumann 
BCs  4>v  =  0  are  used.  The  total  number  of  iterations  needed  for  the  maximum 
change  in  the  solution  in  subsequent  iterations  to  be  reduced  to  <  1CH6  in  this 
example  is  836. 

In  Figure  8  we  show  contour  plots  of  Wulff  shapes  arising  from  solving 


\Jp2  +  q2  +r2 


'2  -L  q2  _|_  y  2 


V/3|r|  -  \Jp2  +  q2 


(p2  +  q2)5/2  -  (-5 p4q  +  10 p2q3  -  q5)  \ 

1  +  V - 2(p2  +  -  =L  (13) 


This  equation  is  equivalent  to  solving 


(1  +  2v/sin(|6»i|  —  0.5tt))(1  +  |  sin(2.5(6»2  +  O.5tt))|)|V0|  =  1. 


The  diffusion  terms  satisfy  a  =  (7/4,  7/4,  2).  The  total  number  of  iterations 
needed  for  the  maximum  change  in  the  solution  in  subsequent  iterations  to  be 
reduced  to  <  10~6  is  429. 

For  these  last  2  examples  the  adaptation  strategy  is  based  on  each  point’s 
I2  distance  to  the  fixed  point,  which  is  not  the  optimal  refinement  strategy 
for  minimizing  the  global  error  in  most  norms.  Thus  there  remains  work  to 
do  in  determining  better  adaptation  strategies.  However,  the  point  of  these 
examples  is  to  show  the  convergence  of  the  fast  sweeping  method  using  the 
LFH  in  a  number  of  sweeps  comparable  to  the  number  found  in  [16]  for  the 
same  examples.  See  [16]  for  more  details  on  these  Wulff  crystal  problems. 


6  Conclusion 

We  have  introduced  a  tree  based  adaptive  method  for  solving  level  set  equations 
which  has  the  advantage  of  being  simplex  free.  Its  implementation  changes  very 
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Figure  7:  Wulff  crystal  with  surface  tension  7  (^pv^yj  =  (1  +  2|  sin 6*1 1) (1  + 
|  sin(1.5(02  +  0.57r))|).  Contours  from  left  to  right  at  </>  =  0.08,0.14,0.2. 


Figure  8:  Wulff  crystal  with  surface  tension  7  =  (1  + 

2y/sin(|0i|  —  0.57t))(1  +  |  sin(2.5(02  +  0.57r))|).  Contours  from  left  to  right  at 
</>  =  0.08,0.14,0.2. 
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little  from  dimension  to  dimension,  allowing  use  by  even  the  most  liyperspatially 
challenged  practitioner.  It  allows  for  well  studied  monotone  numerical  Hamilto¬ 
nians  to  be  used,  avoiding  the  complications  that  often  arise  when  constructing 
monotone  numerical  Hamiltonians  on  unstructured  grids.  The  method  has  been 
applied  to  codimension-m,  1  <  in  <  n,  linear  and  nonlinear,  first  and  second 
order,  static  and  time  dependent  H-J  problems  in  up  to  7cl. 

Future  work  will  include  improving  the  adaptive  procedure  for  static  prob¬ 
lems,  based  on  user  defined  global  errors.  Also,  WENO  type  methods  will  be 
explored  to  increase  the  accuracy.  Although  we  do  not  show  numerical  examples 
here,  the  method  could  be  extended  to  higher  (>  2)  order  nonlinear  PDEs  such 
as  Willmore  flows  [12].  Finally,  it  should  be  noted  that  although  the  methods 
presented  are  done  so  for  H-J  problems,  they  could  be  applied  to  other  applica¬ 
tions  requiring  adaptive  meshes  and  function  interpolation  and  reconstruction. 
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